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We derive conditions under which f{Q) gravity models, whose Lagrangian densities / are written 
in terms of a Gauss-Bonnet term Q, are cosmologically viable. The most crucial condition to be 
' satisfied is d^//d^^ > 0, which is required to ensure the stability of a late-time de-Sitter solution 

' as well as the existence of standard radiation/matter dominated epochs. We present a number of 

, explicit f{Q) models in which a cosmic acceleration is followed by the matter era. We find that the 

■ equation of state of dark energy can cross the phantom divide before reaching the present Universe. 
, ^ ' The viable models have asymptotic behavior d?f/dQ^ — > -1-0 for \Q\ — > oo, in which case a rapid 

, oscillation of perturbations occurs unless such an oscillating degree of freedom is suppressed relative 

. to a homogeneous mode in the early universe. We also introduce an iterative method to avoid 

' numerical instabilities associated with a large mass of the oscillating mode. 

m ; 

I. INTRODUCTION 

The late-time cosmic acceleration can, in principle, originate from a modification of gravity rather than an exotic 
f-H ' source of matter with a negative pressure. Over the past five years, a lot of works on modified gravity have been 
I— —I \ done to identify the origin of dark energy (DE) [l[ . The attractive point in modified gravity models is that they are 
generally more strongly constrained from cosmological observations and local gravity experiments than the models 
' based on the exotic source of matter. 

Presumably the simplest extension to Einstein gravity is the so-called f{R) gravity in which / is an arbitrary 
function of the Ricci scalar R [3]. Even in this simple case it is not generally easy to construct viable f{R) models 
— . , that are consistent with cosmological and local gravity constraints. The main reason for this is that f{R) gravity 

■ gives rise to a strong coupling between DE and a non-relativistic matter in the Einstein frame The models need 
to be carefully designed so that a scalar degree of freedom ("scalaron" is nearly frozen to suppress an effective 
coupling between the scalar field and matter. 

00 \ The conditions for the cosmological viability of f{R) models have been derived in Ref. Q. Among those conditions 
. the requirement, d^f/dR^ > 0, is particularly important to give rise to a saddle matter era followed by a late- 

^ ' time cosmic acceleration. This is also required for the stability of cosmological perturbations @ as well as for the 
• ^ consistency with local gravity experiments [2]- The cosmologically viable f{R) models need to be close to the A-Cold 
^ Dark Matter (ACDM) model in the deep matter era, but the deviation from it becomes important around the late 

^ [ stage of the matter era. Several examples of such viable models were presented in Refs. p, |^. 

- - ' Hu and Sawicki fld\ and Starobinksy [ll[ proposed f{R) models that are consistent with local gravity constraints 
as well as cosmological constraints (see also Refs. [12]). In these models it is possible to find an appreciable deviation 
from the ACDM cosmology. This leaves a number of interesting observational signatures such as the peculiar evolution 
of the DE equation of state [H, [l3| , the modification of the matter power spectrum [ll|, [H, [l^ [3| and the change 
of the convergence spectrum of weak lensing [itI. [l^. This is a welcome feature to distinguish f{R) models from the 
ACDM cosmology in future observations. 

There are other modified gravity DE models that are the generalizations of f{R) gravity. For example Carroll et al. 
[l9t proposed theories with the Lagrangian density f{R, P, Q), where / is an arbitrary function oi R, P = R^i,R>^^ and 
Q = R^vpaR^'^^'^ (here i?^^ and Rf_ii,pa are the Ricci tensor and the Riemann tensor, respectively). These theories are 
plagued b y the appea rance of spurious spin-2 ghosts unless a Gauss-Bonnet (GB) combination, i.e., / = f{R, Q — AP), 
is chosen |2Cll . l2ll . [22l . [23| . Even in this case the graviton itself may still become a ghost in the Friedmann-Lemaitre- 
Robertson- Walker (FLRW) spacetime, unless some no-ghost conditions are verified on the background [23^ . 



o 



'Electronic address; |shinji@rs. kagu.tus.ac.jp] 



2 



The GB curvature invariant Lagrangian, Q = B? — AR^^R^^^ -\-R^^pr,R^^'"^ ^ is a total derivative in the 4-dimensional 
FLRW background. In order to give rise to some contribution of the GB term to the Friedmann equation, we require 
that (i) the GB term couples to a scalar field 0, i.e., F{(f))Q, or (ii) the Lagrangian density / is a function of Q, i.e., 
f{G)- The GB coupling in the case (i) appears in low-energy string effective action ^23] and cosmological solutions 
in such theory have been studied in great details ^]. It was shown by several authors that a late-time cosmic 
acceleration followin g a (sca ling) rn atter era occurs for an exponential coupling (x e^"^ in the presence of a 

scalar-field potential [261. [27l. l28l |29|. [30|. Amendola et al. (3l| studied local gravity constraints in such models and 
showed that the energy contribution coming from the GB term needs to be strongly suppressed for the compatibility 
with solar-system experiments. Thus, in the case (i), it is generally difficult to satisfy local gravity constraints if the 
GB term is responsible for DE. 

In the context oi f(Q) gravity there exists a de-Sitter point that can be used for cosmic acceleration [s^ (see also 
Ref. (33|). It was shown in Ref. that the model with inverse powers of linear combinations of quadratic curvature 
invariants, i.e., f{G) = G" with n < 0, is not cosmologically viable because of the presence of separatrices between 
radiation and DE dominations. In Ref. [s^ it was found that the model f{G) = G^ with n > can be consistent 
with solar-system tests for n ^ 0.074 if the GB term is responsible for DE. Li et al. 36] showed that it is difficult to 
reproduce standard expansion history of the Universe unless f{G) is close to cosmological constant. 

In this paper we present a number of explicit models of f{G) gravity that are cosmologically viable. These models 
mimic the ACDM cosmology in the deep matter era, but the deviation from it becomes important at late times 
on cosmological scales. This situation is similar to the viable f{R) models proposed by Hu and Sawicki 10] and 
Starobinsky pj| . In f{G) gravity, however, the GB term changes its sign during the transition from the matter era 
to the accelerated epoch. We need to take into account this property when we construct viable f{G) models. For 
example, the f{G) model that replaces R in the model of Hu and Sawicki or Starobinsky for G is not cosmologically 
viable. 

There is another difference between f{R) and f{G) theories. The Ricci scalar R vanishes for the vacuum 
Schwarzschild solution, whereas the GB term has a non-vanishing value much larger than Hq around the compact 
objects 'SS] {Hq is the present Hubble parameter). In the presence of matter with a density pm the Ricci scalar R 
is roughly of the order of pm so that one has R/Hq ~ pm/pc, where pc is the present cosmological density. The 
viable f{R) models [13 [HI are designed to have the suppression term (R/Hq) " (n > 0) in addition to the ACDM 
Lagrangian in the region of high density {pm 3> Pc)- In the f{G) gravity this sort of suppression occurs even in 
the vacuum background because of the condition G ^ Hq. Hence the f{G) models might be less constrained by 
local gravity constraints relative to the f{R) models. For the same reason, one could expect that for the interior 
star solutions these modifications of gravity could remain small corrections, provided that singularities of the kind 
f^gg = d^//dC/^ = 0,00 for finite values of G are not encountered. The models discussed in this paper have exactly 
this feature. Note, however, that a detailed study of these issues is needed in order to to further constrain these 
modifications, which we leave for future work. 

We will show that the stability condition for a late-time de-Sitter point is given by f^gg > 0. This can be also derived 
by considering the stability of radiation and matter points. Note that the same condition has been derived in Ref. (36j 
by studying the evolution of cosmological perturbations. Since Li et al. (36j used the metric signature (-1-,— ,— ,— ) 
instead of (—,+,+,+) that we adopt throughout this paper, their stability condition f gg < corresponds to our 
stability condition f^gg > 0. 

For the viable f{G) models we propose in this paper, the second derivative f^gg approaches +0 as \G\ gets larger. 
In this case the perturbations in the Hubble parameter H have a large mass squared proportional to l/(H'^f^gg) 
during radiation and matter eras. This can give rise to rapid oscillations of the Hubble parameter and matter density 
perturbations, as they happen for viable models in f{R) gravity [HI, [l^l- In order to avoid this, the oscillating mode 
needs to be suppressed relative to the homogeneous mode in the early Universe. The numerical instability we typically 
face in radiation and matter epochs is associated with the appearance of this oscillating mode. We shall introduce an 
iterative method to avoid this numerical instability in the high-redshift regime. 



II. COSMOLOGICALLY VIABLE f{g) DARK ENERGY MODELS 

We start with the following action 

S = J d^xV^ [^R + f{G)] + 5„,(.9^,., *,„) , (1) 

where i? is a Ricci scalar, and / is a general differentiable function of G, Sm is a matter action that depends on 
a spacetime metric g^,y and matter fields '^m- We choose units such that = SttGn = 1, where Gjv is a bare 
gravitational constant. It should be pointed out that this action, at the classical level, can be rewritten as an 
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auxiliary scalar field coupled to the Q term, as shown in Ref. '^s'l, following a trick used for the f{R) theory. However 
there is no conformal transformation separating Q from such a field, unlike the f{R) theory in which the conformal 
transformation leads to an Einstein frame action with a canonical scalar field coupled to matter. An important 
quantity in f{Q) gravity is l/f,gg, which plays the role of an effective mass for the theory (as we shall see later). 
The variation of the action ([T]) with respect to g^^ leads to the following equation 

G/iu + 8 [Rf^pua + Rpugcr/i — Rpcr9ufj. — Rfiugap + RpaQup + (ff/ifffcrp ~ Sm^Si^p)] ^''V'^/g + {Gf,g ^ f)9tiu = Tp^ , (2) 

where Gpi, — Rpu — (l/2)i?g^^ is the Einstein-tensor. For the energy momentum tensor Tp^ of a matter fluid we take 
into account the contributions of non-relativistic matter (energy density pm) and radiation (energy density Prad)- In 
a flat FLRW background with the metric ds^ = — d<^ + a(<)^dx^, the 00 component of Eq. ^ gives 



(3) 



where H = a/ a, f g = df/dG, and a dot represents a derivative with respect to cosmic time t. The GB term is given 
by 

g = 24:H^{H'^ + H) . (4) 
The energy densities pm and prad satisfy the continuity equations pm + iH pm — and Prad + 4i/prad = 0, respectively. 

A. Stability of de Sitter point 

Let us first discuss the stability around a de Sitter point present in f{Q) gravity by neglecting the contribution of 
pressure-less matter and radiation. The Hubble parameter, H — Hi, at the de Sitter point satisfies 



3i/?-ai/,e(ai)-/(ei), 



(5) 



where Gi = 2AHf. Note that we used the relations Hi = and Qi = 0. Considering a linear perturbation SHi about 
the de Sitter point, Eq. ([3]) gives 



6Hi = mlf,gg{Hi) HiSg{Hi) - 5g{Hi) 



(6) 



Substituting the relations 5g{Hi) = 24:{AHf5Hi + HfSHi) and Sg{Hi) = 24:H^{6Hi + AHiSHi) into Eq. ®, we 
obtain 



5Hi + 3Hi5Hi 



1 



96i/J>/,ee(iJi) 



HfSHi = . 



This shows that the effective mass squared is [(96i?f /.g;g(i?i)) ^ — 4]iJ^. The solution to Eq. ([7]) is given by 

iHi 



5 Hi — cie 



C2e 



-1± Wl- - 



1 



9 \96Hff^gg{Hi) 



-4 



(7) 



(8) 



where ci and C2 are integration constants. This shows that the de Sitter point is stable under the condition 

< iff (i/i )< 1/384, (9) 

which requires that f^gg{Hi) > 0. Especially when < Hif gg{Hi) < 1/600, the de Sitter point corresponds to a 
stable spiral (damping with oscillations). 

B. Stabilities of radiation and matter points 

We shall also study the stability of matter and radiation fixed points by using a similar method discussed above. Let 
us investigate the case in which the evolution of the scale factor is given hy a oc (p: constant) with the dominance 
of a fluid characterized by an energy density pM (either pm or prad)- Then we have 



3i/2 = gf,g - / - 2AH^f^ggg + PM . 



(10) 
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We consider the first-order perturbations 5h and 5m as follows: 



H = (1 + 5h) , PM = P^i^ (1 + Sm) . 



(11) 



Here the subscript "(b)" represents background values, but in the following we omit it for simplicity. Taking the 
homogeneous perturbations of Eq. (jlOp and using the approximate relation 3H^ ~ pu , we find 



OH 



6 

P 



mH'^f, gggl-p 



l92H^f.gg 



f^gg 



H5 



H 



21(1 -p) 



mH^lgg 



In the limit p — )■ oo without matter perturbations, Eq. ([7]) is recovered. 



1 I ( 1 

f.gg 



The solution to Eq. (fT2)) is described by the sum of the matter-induced mode 5^"^^ and the oscillating mode 



1-p 



(12) 



r(osc) 



[illil 



r(ind) 

-'h 



r(osc) 



(13) 



The matter- induced mode corresponds to a special solution to Eq. induced by the matter perturbation S m ■ The 
oscillating mode is the solution of the equation with Sm = in Eq. ([72]) . For the /{G) models whose deviation from 
the ACDM model is small during radiation and matter eras, f^gg is close to 0. In this case the mass squared 



96H^f,gg ' 



(14) 



is the dominant contribution in front of the term Sh in Eq. (|12p . In order to avoid a violent instability of perturbations 
we require that > 0, giving the condition f^gg > 0. In this case the perturbation 6^^^'^^ oscillates with a frequency 
of the order of M. 

During radiation and matter dominated epochs the GB term evolves a.s G = —24H'^ and G = —12H'^, respectively. 
Since G = 24i?^(a/a) from Eq. Q, the GB term changes its sign at the onset of the late-time acceleration. Hence the 
condition f^gg > needs to be satisfied in the region G < Gi- Since the term 24:H^f^ggG on the r.h.s. of Eq. ([TU]) is of 
the order of H^f,gg, this is suppressed relative to 3H^ under the condition H^f^gg ^ 1. In order for this condition 
to hold in the radiation and matter eras, we require that the term f gg approaches -1-0 with the increase of \G\- 



C. Viable f{g) models 

If we consider a spherically symmetric body (mass Mq and radius r©) with a homogeneous density, it was shown in 
Ref. m that the GB term inside and outside the body is given by ^ = -48{GnMq)^ /r% and G = ■i^iGNM^f /r^ , 
respectively (r is a distance from the center of symmetry). In the vicinity of the Sun or the Earth, is much larger 
than the present cosmological GB term. Go- As we move from the interior to the exterior of the star, the GB term 
crosses from negative to positive. This means that f{G) and its derivatives with respect to G need to be regular for 
both negative and positive values of G whose amplitudes are much larger than Go- 

From the above discussions the viable models need to satisfy the following conditions: 

• (i) f{G) and its derivatives f^g, f,gg,--- are regular. 

• (ii) f,gg > foi' G l£ Gi and f_gg approaches +0 in the limit \G\ — *■ cjo. 

• (iii) < Hlf^gg{Hi) < 1/384 at the dc Sitter point. 
A number of examples for the viable forms of f^gg are 

A 2A A 



where A, n and G* are positive constants. These satisfy the condition f^gg > for all values of G- In the following let 
us study the case (a) with n = 1, the case (b) with n = 2, and the case (c). Integrating f gg with respect to G twice, 
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we obtain the following models 



_arctan(^-j--VG,ln(^l + - 



(A) f{g) = arctan ( ) - - ^G, In ( 1 + ^ ) - aX^G, , (16) 




(B) /(^;) = A-^arctan( ;J ) -aAV^^*, (17) 



(C) f{g) = Av/a7ln cosh (^l^jj - aX^g, , (18) 

where a is a constant. We have dropped the terms proportional to g, since they do not give rise to any contribution 
to the evolution equations. If we demand the condition f[Q — Q) = Q then we have a = 0. These models are also 

consistent with the regularity condition (i). Note that the model f (g) — — A^y^[l — (1 + ^J^/^^) which is a 
generalization of the viable f{R) model proposed by Starobinsky |lll |. is not compatible with the condition f gg > 

for g <Q. 

In the models (A), (B), (C) the term gfg — f on the r.h.s. of Eq. ^ in the regime \g\ ^ C/* can be estimated by (A) 
gf,g - / ~ (A/2) [Hgygl) + 2a] , (B) g/^g - / ~ XVg:il + «), (C) gf,g - / ~ Av^(ln2 + a), respectively As 
long as is of the order of Hq (where Hq is the present Hubble parameter) the term gf^g — / is subdominant relative 
to the term during radiation and matter eras (for aX of the order of unity). Moreover we have H^f.gg ^ 1 for 
1^1 ^ g* in the above three models, which means that the condition \24H^ f ggg\ <C is satisfied in this regime. 
The contribution coming from the f{g) term becomes important when H decreases to the order of H^,. 

Let us discuss the condition (iii) for the model (A) . At the de Sitter point this model satisfies the following relation 

A = : — „ ,„ „, , where yi = —r-n ■ (19) 

2a + ln(l + 2422/8) ' ^ ' 

When a = the r.h.s. has a minimum value Amin = 0.828 at y\ = 0.736. Hence two de Sitter points exist for 
A > 0.828. The stability condition for the de Sitter solution corresponds to 

A < £±^1^1 , (20) 

which gives y\ > 0.736 by using Eq. p^ . Hence one of the de Sitter points that exists in the region yi > 0.736 is 
stable. 

In the case where a — we have numerically found instabilities of cosmological solutions around the region g = 
(which occurs during the transition from the matter era to the accelerated era). One can estimate the stability of 
solutions by setting p = 1 in Eq. . This gives the stability condition 

< Hi f,gg{H2) <1/38A, (21) 

where H2 is the Hubble parameter at g = 0. For the model (A) this translates into 

A < -^—rr , (22) 

384y6 ' ^ ' 

where y2 = H2/gV'^. Here 2/2 is slightly larger than j/i. Since yi > 0.736 for a — 0, the condition (j22p requires that 
A ^ 1. However this is not compatible with the condition A > 0.828 for the existence of de Sitter solutions. As we will 
see in the next section, it is generally difficult to have a natural transition around g = for a — 0. This anticipates 
that the cosmological constant term may be needed in general for the cosmological viability of the f{g) models. If 
this is the case, the need of such a constant term makes the f{g) models less attractive from a theoretical point of 
view. 

If a 7^ it is possible to make A much smaller than 1 provided that a 3> 1 [see Eq. ([T5| ]. When a — 10^ and 
2/1 = 0.1, for example, we have A = 3.0 x 10"'* and hence the condition is satisfied even for 2/2 = 1. In the next 
section we shall show that viable cosmological trajectories can be realized for a larger than the order of unity. 

The stabilities of cosmological solutions for the models (B) and (C) are similar to those for the model (A) discussed 
above. The difference is that the models (B) and (C) have larger powers of g"^ /g1 in f gg compared to the model (A). 
Since the mass grows rapidly toward the past in such models, it is more difficult to start solving the equations 
numerically from the high-rcdshift regime unless we use an iterative method we discuss later. 
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D. Oscillating modes in the early Universe 



At the end of this section we discuss the evolution of the homogeneous perturbation 6h during radiation and matter 
eras for the viable f(G) models presented above. For the matter-induced mode, the mass term M'^S^ balances with 
the source term on the r.h.s. of Eq. (|12p. giving 



c(ind) 



Sm/2 . 



(23) 



Hence the matter- induced mode grows in proportion to 5m ■ In the matter era 5^"'^^ oc t^/"^ in the regime where the 
model is close to ACDM model [13 . 

Let us consider the model (B) for the evolution of the oscillating mode in the regime \G\ ^ G*- In this case we have 
9QH'^{f^ggg / f^gg){l — p)/p^ — 16/p and hP ~ fi^t^^^ , where ^ is a constant. Hence the oscillating mode satisfies 



•• (osc) 3p+ 10 • (osc) jJi r-(osc) 

iH H ; Oh + T^^h 



0. 



t " ti 

The solution to this equation can be written as the combination of Bessel differential functions: 

l+5/p 



j(osc) _ / U 

" ' t 



Je, {z) + 



^{2} 
J-Si (Zi) 



J-eAz) 



(24) 



(25) 



where 9i ~ (5 +p)/{5p), z = fj,/{5t^), z\ = z(ti), and ti is the initial time at which two modes in the square bracket 
of Eq. 
As t 



25)) have amplitudes ^^"^ ^^\- 



0, the solution reduces to 



5f'^ ~ ^i-(10-3P)/(2p)<!Ci 



(l + 20i) 



t 



C2 cos 



(26) 



where A 



-1/4 



•^/TO/tt, and Ci and C2 are constants. During the matter era [p — 2/3), the amplitude of 5)^ 



(osc) 



is proportional io t ^ so that the oscillating mode tends to be negligible relative to the matter-induced mode with 
time. However, as we go back to the past, the amplitude grows with larger frequency of oscillations. Since one has 
Am p|"(5^ ^''^] cx for p = 1/2, this property also persists during the radiation-dominated epoch. As we see in 

Sec. lHII the large mass term M tends to lead numerically instabilities associated with violent oscillations of 5h, unless 
the oscillating mode is strongly suppressed relative to the matter-induced mode. For the model (A) the growth of 
the mass squared is not so strong (M^ oc t^^) relative to the model (B), but still it is difficult to solve equations 
numerically from the high-redshift regime. This property is even severe for the model (C) because of the rapid increase 
of the mass squared: cx i'' exp(c/i'') (c is a positive constant). 



III. COSMOLOGICAL DYNAMICS 



In this section we discuss cosmological dynamics for the viable models presented in the previous section. In 
subsection IIII Al we first integrate the background equations for the model (A) directly in the low-redshift regime. It 
is more difficult to solve the equations numerically as we start integrating from higher redshifts. This is related to the 
appearance of the oscillating mode with a large frequency. In subsection IIIIBI we shall propose an iterative method 
to get approximate solutions in such a situation. 



A. Low-redshift cosmological solutions 

In order to discuss cosmological solutions in the low-redshift regime, it is convenient to introduce the following 
dimensionless quantities 

= — :^ — O = O = /^rad ,„„s 
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where = G 



1/4 



We then obtain the following equations of motion 



-4a; - 4a; 



1 



Qf,g - f 



3(1 ^rad) 



y =xy , 

n'^ = -(3 + 2a;)r!„, 
^Ld - -(4 + 2x)aad 



(28) 

(29) 
(30) 
(31) 



where a prime represents a derivative with respect to = ln(a). The quantities H^f^gg and {QJ^g — f)/H'^ can be 
expressed by x and y once the model is specified. The energy fraction of DE is given by Que = 1 — ~ i^rad- We 
also define the effective equation of state 



2H 2 

which changes from to —1 from the matter era to the final de Sitter epoch for viable f{Q) models. 

2.0 

1.0 

0.0 
-1.0 
-2.0 
-3.0 
-4.0 



(32) 
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FIG. 1: The evolution of Wcff, SIde and Qm versus the scale factor a for the model (A) with parameters a = 10.0 and 
A — 7.5 X 10~^. The initial conditions are chosen to be a; = —1.502, y = 20.0, fim = 0.9959 and firad = 0.004. These results 
are obtained by integrating Eqs. (|28|l - (|3ip forward using the Isode stiff integrator. It is clear that the matter era is followed by 
the accelerated epoch with the oscillation of Wcs around —1. 

In Fig.[T]we plot the evolution of Wcs, ^be and for the model (A) with parameters a — 10.0 and A — 7.5 x 10~^. 
The present epoch corresponds to the scale factor a — 1 with JIde = 0.72 and flm = 0.28. From Eq. (fTg]) there exists 
a de Sitter point at yi — 0.518 that satisfies the stability condition pp)) . In fact Fig. [2] shows that the quantity 
y — H/Hir approaches this value in the asymptotic future. 

In Fig.[T]the effective equation of state Wcff oscillates around —1 as the system enters the epoch of cosmic acceleration, 
which implies that the de Sitter solution is a stable spiral. It is interesting to note that w^s drops down to a value less 
than —4 around the present epoch. We also find in Fig. [5] that the GB term switches its sign during the transition 
from the matter era to the accelerated epoch, (which corresponds to passing through the minus infinity in logarithmic 
scale). 

We have also tried numerical integrations by changing the model parameters a and A. For the values of a smaller 
than unity it is not easy to to get plausible cosmological evolution. This is associated with the fact that the stability 
condition (|22p is difficult to be satisfied around = for smaller a. In Fig. [3] we illustrate the variation of Wcff, f^DE 
and for the model (A) with parameters a = and A = 1. While a stable de Sitter point exists at yi — 1.075, 
VLm becomes larger than the order of unity during the transition from the matter era to the accelerated epoch. This 
reflects the instability of the solutions around Q = Q. For a smaller than the order of unity, the solutions exhibit 
unusual behavior similar to that in Fig. [3] or they do not reach the de Sitter attractor. 
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FIG. 2: The evolution of the quantities H/Ht,, \Q\ and H^f^gg for the same model parameters as given in Fig. [T] The GB 
term changes its sign from negative to positive during the transition from the matter era to the accelerated epoch. The term 
H^f,gg becomes much smaller than 1 as we go back to the past. 




FIG. 3: The evolution of Wcft , f^DE and ilm versus the scale factor a for the model (A) with parameters a — (i.e., f{Q = 0) = 0) 
and A = 1. The initial conditions are chosen to be a; = —1.502, y = 20.0, Sim = 0.9959 and f^iad = 0.004. From the matter era 
to the accelerated epoch the solution shows an unusual transition where Urn exceeds the order of unity. 



For the model (A) it is difficult to integrate the equations from the redshift larger than 50 because the term 
l/{96H^f,gg) in Eq. (fT2|) gets very large as we go back to the past. The decrease of the term H^f^gg for larger z is 
in fact confirmed in Fig. [21 On one hand, this is a good property, as the spin-2 no-ghost conditions are then satisfied 
[2^, in addition to the fact that scalar perturbations remain stable [3^. However, this large mass leads to a rapid 
oscillation of the perturbations Sh, which is hard to be dealt with numerically. The difficulty of numerical integration 
is even more severe for the models (B) and (C). When the equations are integrated forward, we need to choose initial 

conditions carefully so that the oscillating mode J^'^'^'' is suppressed relative to the matter- induced mode (J^"*^"*. This 
property is similar to viable f{R) gravity models in which the system is unstable in the high-redshift regime because 
of the appearance of oscillating mode 11, 13]. 
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B. High-redshift approximate solutions 



In the regime of high-redshifts one can use an iterative method (known as the "fixed-point" method) to find 
cosmological solutions approximately. We define and G to be = H^/Hq and G = G/Hq, where the subscript 
"0" represents present values (with gq = 1). The models (A), (B) and (C) can be written in the form 

/(a)-/(a)i/o'-Ai/2, (33) 

where A — aX^/G^ / Hq and f{G) is a function of G- The modified Friedmann equation ([3]) reduces to 

H'-Hl = l if.gQ - I) (34) 



where 



Note that Hp^ represents the Hubble parameter in the ACDM model. In the following we omit the tilde for simplicity. 
In Eq. ([51)1 there are derivatives of H in terms of N up to the second-order. Then we write Eq. ([51)1 in the form 



H^-HI = c(h\H^\h^") , (36) 

where C = [f^gG - /)/3 - Si/'' (d/^g/dTV). At high redshifts (a < 0.01) the models (A), (B) and (C) are close to 
the ACDM model, i.e., ~ H\. We shall introduce an iterative method to derive approximate solutions in such a 
regime. 

As a starting guess we set the solution to be H"^!^^ = H\. The first iteration is then 

^fi) = + C(o) , (37) 

where C(o) = C^i^H'^^y H'^q^' , H'^^^"^ . This first iterative solution was used as an approximate solution for inverse 
curvature gravity in the paper of Mena et al. {3^ , while the authors did not pursue the idea of iterating the process 
again. We shall iterate the process in order to get better approximate solutions. The second iteration is 

^(2) ^ H\ + C(i) , (38) 

where C(i) = C{Hf,y Hf,^' , Hf,^") . 

If the starting guess was in the basin of a fixed point, H'^^^ will converge to the solution of the equation after the 
z-th iteration. For the convergence we need the following condition 

which means that each correction decreases for larger i. The following relation is also required to be satisfied: 

Once the solution begins to converge, one can stop the iteration up to the required/ available level of precision. 

At very high redshifts (say, the epoch of nucleosynthesis), the above method is presumably the only one that 
provides approximate cosmological solutions. We have checked that, for N > —4 (i.e., for the redshift z < 50), this 
iterative method and the direct-forward-integration give the same results. To be more precise, the iteration is used 
in order to find the values of and G 'd-t N = —4. Then these values are adopted as initial conditions for the 
direct-forward-integrator. We integrate the equations for —4 < N < —3 (where the iterative method still works well) 
and compare H"^ as well as ^ at = — 3 derived by two methods. We find that these provide identical results with 
the precision of the order of 10~^. This shows that the solution derived by direct integration remains close to the 
iterative one at least for all values of N at which the initial guess for the iterative method is in the basin of the fixed 
point. 
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FIG. 4: The plot of the absolute errors log^gdHf 



a I) (left) and log^o 



Hf-Hf+C, 



^] (right) 



versus A*' for the model 



(A) with i = 0, 1, • ■ ■ ,6. The model parameters a and A are the same as those in Fig. [T] The iterative method provides the 
solutions with high accuracy in the regime A'^ ^ —4. 
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model (A) with the same model parameters as used in[T] 
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Ci\) as well as the relative error logj^ 



In Fig. [3] we plot the absolute error logj^gd-^i 

model (A) with a ~ 10.0 and A — 7.5 x 10~^ (i.e., the same model parameters as used in Fig. [Ij 
have carried out the iteration for 6 times. The absolute error logj^Q 
Hf — H\ is really close to C, 



\Hf-Hi-C,\ 



for the 



\m -Hi -a 



_\H'f-Hl+Ci 

Note that we 
is not sufficient to confirm that 
However the smallness of the relative error in the left panel of Fig. [5] confirms that 
the solution derived by the iterative method is very accurate. While this approximation tends to be worse for lower 
redshifts, the direct integration is well suited for N ^ —4 as we presented in the subsection IIII Al 

The left panel of Fig. O shows that, for N < —4, the iterative solution is very similar to the ACDM solution 
characterized by Hi. Hence the Universe passes through the radiation-dominated epoch to the matter-dominated 
one as in the ACDM cosmology. From the right panel of Fig. [5] we find that the quantity fggH^ is very much smaller 
than unity for A'^ ^ —4, which leads to an extremely large frequency M for the perturbation Sh- This is the main 
reason why we need to use the iterative method rather than the direct integration to avoid numerical instabilities in 
the high-redshift regime. 



IV. CONCLUSIONS 



In this paper we have constructed viable f{G) gravity models that are cosmologically viable. In order to have a 
stable de Sitter attractor the condition ^ needs to be satisfied. For the stability of radiation and matter fixed points 
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the mass squared given in Eq. (|14p is required to be positive. These results show that the quantity f^gg must 
be positive to obtain viable cosmological evolution. Since the GB term changes its sign during the transition from 
the matter era to the epoch of cosmic acceleration, we need to construct models in which neither the violation of the 
stability conditions nor the divergence of some terms occurs in the past expansion history of the Universe. 

A number of explicit f{G) models satisfying the above requirements are given in Eqs. p6l) - (fT8)) . Even if these 
models do not require an exotic source of matter responsible for the cosmic acceleration, it should be pointed out that 
in the examples presented here a cosmological constant term is still required. These models come from the integration 
of viable forms of f,gg presented in Eq. (jlSp . In the regime where \G\ is much larger than t/* (which is the same order 
as the present value Go), the model mimics the ACDM cosmology. The deviation from the ACDM cosmology tends to 
be important a,s \G\ approaches the order of G*- Since the mass squared becomes very much larger than as we 
go back to the past, this leads to rapid oscillations of the Hubble parameter unless initial conditions are chosen such 
that the oscillating mode is suppressed relative to the matter-induced mode. The direct integration of Eqs. (|28p - ([3T|) 
is prone to numerical instabilities in the high-redshift regime because of the very heavy mass M . 

We have adopted an iterative method to derive approximate solutions in the high-redshift regime. The results of 
Figs [4] and [5] show that the iterative method gives rise to the cosmic expansion history that is very close to the ACDM 
model for z > 50. We have used these results as initial conditions at 2; ~ 50 for the direct forward- integration in the 
low-redshift region. We have found that the effective equation of state WcS enters the phantom region {wcS < —1) 
before reaching the de Sitter attractor with Wcs = —1 (see Fig. [Ij. For the models with f(G = 0) = the solutions 
typically exhibit unusual behavior where grows larger than 1 during the transition from the matter era to the 
accelerated era (see Fig. [3]). This is associated with the fact that an instability around = is present for small 
a (< 1) for the models P^ - (|18p . while this is not the case for a larger than the order of unity. 

Since the f{G) models we have proposed mimic the ACDM model in the high-curvature regime whose energy density 
is much larger than the present cosmological one, it should be possible for them to satisfy local gravity constraints. 
We leave detailed analysis for the compatibility of our models with local gravity experiments for future work. 
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